Libraries

library(survival)
library(KMsurv)
library(ggplot2)
library(ggpubr)
library(survminer)
library(plotly)
library(muhaz)
library(ggthemes)

Loading in Data

data(burn)

burn1 <- burn
burn1 <-
  data.frame(burn1, Treatment = factor(burn1$Z1, labels = c("Routine", "Cleansing")))
burn1 <-
  data.frame(burn1, Gender = factor(burn1$Z2, labels = c("Male", "Female")))
burn1 <-
  data.frame(burn1, Race = factor(burn1$Z3, labels = c("Nonwhite", "White")))
burn1 <- data.frame(burn1, PercentBurned = burn1$Z4)
burn1 <-
  data.frame(burn1, SiteHead = factor(burn1$Z5, labels = c("NotBurned", "Burned")))
burn1 <-
  data.frame(burn1, SiteButtock = factor(burn1$Z6, labels = c("NotBurned", "Burned")))
burn1 <-
  data.frame(burn1, SiteTrunk = factor(burn1$Z7, labels = c("NotBurned", "Burned")))
burn1 <-
  data.frame(burn1, SiteUpperLeg = factor(burn1$Z8, labels = c("NotBurned", "Burned")))
burn1 <-
  data.frame(burn1, SiteLowerLeg = factor(burn1$Z9, labels = c("NotBurned", "Burned")))
burn1 <-
  data.frame(burn1, SiteRespTract = factor(burn1$Z10, labels = c("NotBurned", "Burned")))
burn1 <-
  data.frame(burn1, BurnType = factor(burn1$Z11, labels = c("Chemical", "Scald", "Electric", "Flame")))

burn1.surv <- with(burn1, Surv(T3, D3))
# plot(survfit(burn1.surv ~ Treatment, data = burn1),
#      col = 1:2,
#      lwd = 2)
# title("Time to Infection for Routine Care and Total Body Cleansing")
# legend(
#   "topright",
#   c("Routine Care", "Total Body Cleansing"),
#   col = 1:2,
#   lwd = 2
# )

#Lets observe the data 
burn1
print(survdiff(burn1.surv ~ Treatment, data = burn1))
## Call:
## survdiff(formula = burn1.surv ~ Treatment, data = burn1)
## 
##                      N Observed Expected (O-E)^2/E (O-E)^2/V
## Treatment=Routine   70       28     21.4      2.07      3.79
## Treatment=Cleansing 84       20     26.6      1.66      3.79
## 
##  Chisq= 3.8  on 1 degrees of freedom, p= 0.05

Question 1

# Plot the Kaplan Meier Curves 
KMcurves <- survfit(burn1.surv ~ Treatment, data = burn1)

KMplot <-
  ggsurvplot(KMcurves, ggtheme = theme_fivethirtyeight()) + labs(title = "KM Curves for Time to Infection for Routine Care \n and Total Body Cleansing")
ggplotly(KMplot[[1]])

Question 2

cumHazPlot <-
  ggsurvplot(
    KMcurves,
    fun = "cumhaz",
    conf.int = FALSE,
    palette = c("#5d8aa8", "#e32636"),
    ggtheme =theme_fivethirtyeight()
  ) + ggtitle("Cumulative Hazard for Treatment Type")

#ggplotly(cumHazPlot[[1]])

#cumHazPlot
ggplotly(cumHazPlot$plot) # + geom_abline())
## Warning: plotly.js does not (yet) support horizontal legend items 
## You can track progress here: 
## https://github.com/plotly/plotly.js/issues/53
#Cloglog plot
cLogLogPlot <-
  ggsurvplot(
    KMcurves,
    fun = "cloglog",
    palette = c("#5d8aa8", "#e32636"),
    ggtheme = theme_fivethirtyeight()
  ) + ggtitle("Complimentary Log-Log for Treatment") 
ggplotly(cLogLogPlot[[1]])
## Warning: plotly.js does not (yet) support horizontal legend items 
## You can track progress here: 
## https://github.com/plotly/plotly.js/issues/53

Question 3

#Building with time independent objects
survTimeIndep <- Surv(burn1$T3, burn1$D3)
#Build a cox model with just treatment
burn.cox.indep.T <- coxph(survTimeIndep ~ Treatment, data = burn1)
summary(burn.cox.indep.T)
## Call:
## coxph(formula = survTimeIndep ~ Treatment, data = burn1)
## 
##   n= 154, number of events= 48 
## 
##                       coef exp(coef) se(coef)      z Pr(>|z|)  
## TreatmentCleansing -0.5614    0.5704   0.2934 -1.914   0.0557 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing    0.5704      1.753     0.321     1.014
## 
## Concordance= 0.566  (se = 0.039 )
## Likelihood ratio test= 3.73  on 1 df,   p=0.05
## Wald test            = 3.66  on 1 df,   p=0.06
## Score (logrank) test = 3.76  on 1 df,   p=0.05
#Building with treatement and burnpercentage
burn.cox.indep.TPB <- coxph(survTimeIndep ~ Treatment + PercentBurned, data = burn1)
summary(burn.cox.indep.TPB)
## Call:
## coxph(formula = survTimeIndep ~ Treatment + PercentBurned, data = burn1)
## 
##   n= 154, number of events= 48 
## 
##                         coef exp(coef)  se(coef)      z Pr(>|z|)  
## TreatmentCleansing -0.524764  0.591695  0.295769 -1.774    0.076 .
## PercentBurned       0.007248  1.007275  0.007145  1.015    0.310  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing    0.5917     1.6901    0.3314     1.056
## PercentBurned         1.0073     0.9928    0.9933     1.021
## 
## Concordance= 0.585  (se = 0.046 )
## Likelihood ratio test= 4.7  on 2 df,   p=0.1
## Wald test            = 4.72  on 2 df,   p=0.09
## Score (logrank) test = 4.82  on 2 df,   p=0.09
#Building with treatement, burnpercentage, and burn Type 
burn.cox.indep.TPBbt <- coxph(survTimeIndep ~ Treatment + PercentBurned + BurnType, data = burn1)
summary(burn.cox.indep.TPBbt)
## Call:
## coxph(formula = survTimeIndep ~ Treatment + PercentBurned + BurnType, 
##     data = burn1)
## 
##   n= 154, number of events= 48 
## 
##                         coef exp(coef)  se(coef)      z Pr(>|z|)  
## TreatmentCleansing -0.537948  0.583945  0.300822 -1.788   0.0737 .
## PercentBurned       0.008986  1.009027  0.007301  1.231   0.2184  
## BurnTypeScald       1.057604  2.879465  1.085913  0.974   0.3301  
## BurnTypeElectric    2.265027  9.631385  1.083574  2.090   0.0366 *
## BurnTypeFlame       0.890441  2.436204  1.019409  0.873   0.3824  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing    0.5839     1.7125    0.3238     1.053
## PercentBurned         1.0090     0.9911    0.9947     1.024
## BurnTypeScald         2.8795     0.3473    0.3428    24.190
## BurnTypeElectric      9.6314     0.1038    1.1517    80.543
## BurnTypeFlame         2.4362     0.4105    0.3304    17.965
## 
## Concordance= 0.654  (se = 0.044 )
## Likelihood ratio test= 12.95  on 5 df,   p=0.02
## Wald test            = 14.22  on 5 df,   p=0.01
## Score (logrank) test = 15.56  on 5 df,   p=0.008
#Building with treatement, burnpercentage, and burn Type + Respartory Tract
burn.cox.indep.TPBbtResp <- coxph(survTimeIndep ~ Treatment + PercentBurned + BurnType + SiteRespTract, data = burn1)
summary(burn.cox.indep.TPBbtResp)
## Call:
## coxph(formula = survTimeIndep ~ Treatment + PercentBurned + BurnType + 
##     SiteRespTract, data = burn1)
## 
##   n= 154, number of events= 48 
## 
##                          coef exp(coef)  se(coef)      z Pr(>|z|)  
## TreatmentCleansing  -0.536426  0.584835  0.300473 -1.785   0.0742 .
## PercentBurned        0.008111  1.008144  0.007588  1.069   0.2851  
## BurnTypeScald        1.083037  2.953636  1.087820  0.996   0.3194  
## BurnTypeElectric     2.278658  9.763566  1.084165  2.102   0.0356 *
## BurnTypeFlame        0.851236  2.342540  1.024638  0.831   0.4061  
## SiteRespTractBurned  0.140348  1.150675  0.351875  0.399   0.6900  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                     exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing     0.5848     1.7099    0.3245     1.054
## PercentBurned          1.0081     0.9919    0.9933     1.023
## BurnTypeScald          2.9536     0.3386    0.3503    24.906
## BurnTypeElectric       9.7636     0.1024    1.1662    81.743
## BurnTypeFlame          2.3425     0.4269    0.3144    17.453
## SiteRespTractBurned    1.1507     0.8691    0.5773     2.293
## 
## Concordance= 0.649  (se = 0.044 )
## Likelihood ratio test= 13.11  on 6 df,   p=0.04
## Wald test            = 14.27  on 6 df,   p=0.03
## Score (logrank) test = 15.66  on 6 df,   p=0.02
#Building with treatement, burnpercentage, and burn Type + resp Tract and gender
burn.cox.indep.TPBbtRespG <- coxph(survTimeIndep ~ Treatment + PercentBurned + BurnType + SiteRespTract + Gender, data = burn1)
summary(burn.cox.indep.TPBbtRespG)
## Call:
## coxph(formula = survTimeIndep ~ Treatment + PercentBurned + BurnType + 
##     SiteRespTract + Gender, data = burn1)
## 
##   n= 154, number of events= 48 
## 
##                          coef exp(coef)  se(coef)      z Pr(>|z|)  
## TreatmentCleansing  -0.590227  0.554202  0.302088 -1.954   0.0507 .
## PercentBurned        0.007192  1.007218  0.007599  0.946   0.3439  
## BurnTypeScald        1.078282  2.939625  1.088765  0.990   0.3220  
## BurnTypeElectric     2.157086  8.645907  1.087701  1.983   0.0473 *
## BurnTypeFlame        0.843227  2.323853  1.024277  0.823   0.4104  
## SiteRespTractBurned  0.169939  1.185232  0.356008  0.477   0.6331  
## GenderFemale        -0.549411  0.577290  0.397753 -1.381   0.1672  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                     exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing     0.5542     1.8044    0.3066     1.002
## PercentBurned          1.0072     0.9928    0.9923     1.022
## BurnTypeScald          2.9396     0.3402    0.3480    24.834
## BurnTypeElectric       8.6459     0.1157    1.0256    72.889
## BurnTypeFlame          2.3239     0.4303    0.3121    17.301
## SiteRespTractBurned    1.1852     0.8437    0.5899     2.381
## GenderFemale           0.5773     1.7322    0.2647     1.259
## 
## Concordance= 0.682  (se = 0.043 )
## Likelihood ratio test= 15.22  on 7 df,   p=0.03
## Wald test            = 16  on 7 df,   p=0.03
## Score (logrank) test = 17.39  on 7 df,   p=0.02